In [1]:
%matplotlib
In [2]:
#%pylab inline
import numpy as np
from matplotlib import pyplot as plt
import os,sys
import scipy as sp
In [3]:
import ewc
import wrm_vangenuchten
import permafrost_model_explicit_fpd
import permafrost_model_explicit_fpd_smoothed
import capillary_pressure
In [4]:
reload(permafrost_model_explicit_fpd_smoothed)
reload(permafrost_model_explicit_fpd)
Out[4]:
In [5]:
colors = dict(A='b', A2='c', B='r', B2='m', C='g', C2='lime')
In [6]:
wrms = dict()
colors = dict()
# mineral soil
mA = 0.2481
aA = 0.00033
srA = 0.199
wrms['mineral soil'] = wrm_vangenuchten.VanGenuchten(m=mA, alpha=aA, sr=srA, smoothing_interval=0.1)
colors['mineral soil'] = 'k'
# peat
mB = 0.3055
aB = 0.00095
srB = 0.34
wrms['peat, Atchley GMD'] = wrm_vangenuchten.VanGenuchten(m=mB, alpha=aB, sr=srB, smoothing_interval=0.05)
colors['peat, Atchley GMD'] = 'r'
# Moon and Graham peat
mD = 0.19
aD = 0.00051
srD = 0.
wrms['peat, Moon and Graham'] = wrm_vangenuchten.VanGenuchten(m=mD, alpha=aD, sr=srD, smoothing_interval=0.05)
colors['peat, Moon and Graham'] = 'm'
# Letts hemic peat
mE = 0.41
aE = 0.000204
srE = 0.17
wrms['peat, Letts hemic'] = wrm_vangenuchten.VanGenuchten(m=mE, alpha=aE, sr=srE, smoothing_interval=0.05)
colors['peat, Letts hemic'] = 'b'
mF = 0.47
aF = 0.000817
srF = .04
wrms['peat, Letts fibric'] = wrm_vangenuchten.VanGenuchten(m=mF, alpha=aF, sr=srF, smoothing_interval=0.05)
colors['peat, Letts fibric'] = 'c'
# moss
mC = 0.2754
aC = 0.0023
srC = 0.05
wrms['moss'] = wrm_vangenuchten.VanGenuchten(m=mC, alpha=aC, sr=srC, smoothing_interval=0.01)
colors['moss'] = 'g'
In [7]:
saturations = np.linspace(0,1,1000)
spcs = dict()
for key, wrm in wrms.iteritems():
sA = saturations[np.where(saturations > wrm._sr)[0]]
pcsA = np.array([wrm.capillaryPressure(s) for s in sA])
spcs[key] = sA,pcsA
In [9]:
# plot the WRMs
plt.figure(figsize=(16,10))
for key, (s,pc) in spcs.iteritems():
#if not key.startswith("moss") and not key.startswith("mineral"):
plt.semilogy(s,pc,color=colors[key], label=key)
plt.legend(loc="upper right")
plt.show()
In [22]:
Ts1 = np.arange(273.15-20, 273.15-5, 1.)
Ts2 = np.arange(273.15-5, 273.15-1, .1)
Ts3 = np.arange(273.15-1, 273.15+1, .001)
Ts4 = np.arange(273.15+1, 273.15+5, .1)
Ts5 = np.arange(273.15+5, 273.15+50, 1.)
Ts = np.concatenate((Ts1, Ts2, Ts3, Ts4, Ts5))
def ewc_func(wrm, poro=0.5, sl0=0.9, left=-1e8, right=1e6, coef=1,
pm_mod=permafrost_model_explicit_fpd, pc_ice=capillary_pressure.pc_ice):
# create the permafrost model
my_wrm = wrm_vangenuchten.VanGenuchten(m=wrm._m, alpha=wrm._alpha*coef, sr=wrm._sr, smoothing_interval=(1-wrm._s0))
pm = pm_mod.PermafrostModel(my_wrm, pc_ice=pc_ice)
if sl0 < 1.:
# determine the water content at sl0, T=273.15
wc0 = ewc.water_content(273.15, 101325. - my_wrm.capillaryPressure(sl0), poro, pm)
# for each temperature, determine the pressure that matches that wc
good_Ts = []
good_ps = []
for T in Ts:
def wcRoot(p):
return ewc.water_content(T,p,poro,pm) - wc0
try:
myp,r = sp.optimize.brentq(wcRoot, left, right, full_output=True)
except ValueError:
print "Bracket failure at T = ", T, "wc = ", wc0
print " wc(", left, ") = ", wcRoot(left)+wc0
print " wc(", right,") = ", wcRoot(right)+wc0
else:
if r.converged:
good_Ts.append(T)
good_ps.append(myp)
e_wc = np.array([ewc.ewc(T,p,poro,pm) for T,p in zip(good_Ts,good_ps)])
sats = np.array([pm.saturations_Tp(T,p) for T,p in zip(good_Ts, good_ps)])
else:
good_Ts = Ts
good_ps = 101325.*np.ones(Ts.shape,'d')
e_wc = np.array([ewc.ewc(T,101325.,poro,pm) for T in Ts])
sats = np.array([pm.saturations_Tp(T,101315.) for T in Ts])
result = dict()
result['T'] = np.array(good_Ts)
result['p'] = np.array(good_ps)
result['e'] = e_wc[:,0]
result['wc'] = e_wc[:,1]
result['porosity'] = poro
result['si'] = sats[:,2]
result['sl'] = sats[:,1]
result['sg'] = sats[:,0]
result['kr'] = np.array([my_wrm.k_relative(sl) for sl in sats[:,1]])
print sats[:,1]
print result['kr']
return result
def surf_func(transition=1.):
def uf(temp):
adj = temp - 273.15
if adj > 0.:
return 1.
elif adj < -transition:
return 0.
else:
return (np.sin(np.pi/2. * (adj + transition/2.)/(transition/2.)) + 1.)/2.
result = dict()
result['T'] = Ts.copy()
result['p'] = 101325*np.ones(Ts.shape,'d')
result['sg'] = np.zeros(Ts.shape,'d')
result['sl'] = np.array([uf(T) for T in Ts])
result['si'] = 1.-result['sl']
result['e'] = None
return result
In [23]:
def getAxs(fig=None):
if fig is None:
fig = plt.figure(figsize=(14,4))
axs = []
axs.append(fig.add_subplot(141))
axs.append(fig.add_subplot(142))
axs.append(fig.add_subplot(143))
axs.append(fig.add_subplot(144))
axs[0].set_xlabel("Temperature [K]")
axs[1].set_xlabel("Temperature [K]")
axs[2].set_xlabel("Temperature [K]")
axs[0].set_ylabel("pressure [Pa]")
axs[1].set_ylabel("energy [J/m^3]")
axs[2].set_ylabel("ice/liq saturation [-]")
axs[3].set_ylabel("relative perm [-]")
return axs
def plot(res, color, name, axs):
axs[0].plot(res['T'],res['p'],'-x', color=color, label=name)
if res['e'] is not None:
axs[1].plot(res['T'],res['e'],'-x', color=color)
#axs[2].plot(res['T'],res['si'],'-x',color=color)
axs[2].plot(res['T'],res['sl'],'--x',color=color)
axs[3].semilogy(res['sl'],res['kr'],'--x',color=color)
def finishPlot(axs):
axs[0].legend(loc="lower right")
In [30]:
print wrmA._spline.__dict__
print wrmA.d_k_relative(0.89999999)
print wrmA.d_k_relative(0.90000000001)
print wrmA._spline.dy1
resultA = ewc_func(wrmA, poro=0.7, sl0=1, pm_mod=permafrost_model_explicit_fpd)
resultA2 = ewc_func(wrmA2, poro=0.7, sl0=1, pm_mod=permafrost_model_explicit_fpd, pc_ice=capillary_pressure.pc_ice_regularized(0.2))
resultA3 = ewc_func(wrmA, poro=0.7, sl0=1, pm_mod=permafrost_model_explicit_fpd, pc_ice=capillary_pressure.pc_ice_regularized(0.2))
#resultA2 = ewc_func(wrmA, poro=0.5, sl0=0.9, coef=.01)
In [ ]:
resultB = ewc_func(wrmB, poro=0.85, sl0=1, pm_mod=permafrost_model_explicit_fpd)
resultB3 = ewc_func(wrmB, poro=0.85, sl0=1, pm_mod=permafrost_model_explicit_fpd, pc_ice=capillary_pressure.pc_ice_regularized(0.2))
#resultB2 = ewc_func(wrmB,0.85, sl0=0.9, right=1e10, coef=10)
In [ ]:
resultC = ewc_func(wrmC,0.9,0.9)
resultC2 = ewc_func(wrmC,0.9,0.9, pm_mod=permafrost_model_explicit_fpd_smoothed)
#resultC2 = ewc_func(wrmC,0.9,0.5, coef=10)
In [25]:
%matplotlib
axs = getAxs()
plot(resultA, 'b', 'soil, fpd', axs)
plot(resultA2, 'c', 'soil, fpd smoothed', axs)
plot(resultA3, 'g', 'soil, fpd smoothed', axs)
finishPlot(axs)
plt.tight_layout()
plt.show()
In [ ]:
axs = getAxs()
plot(resultB, 'r', 'peat, fpd', axs)
plot(resultB3, 'm', 'peat, fpd smoothed', axs)
finishPlot(axs)
plt.tight_layout()
In [ ]:
axs = getAxs()
plot(resultC, 'g', 'moss, fpd', axs)
plot(resultC2, 'lime', 'moss, fpd smoothed', axs)
finishPlot(axs)
plt.tight_layout()
plt.show()
In [ ]:
axs = getAxs()
plot(resultA, 'b', 'soil, s=0.9', axs)
plot(resultA2, 'c', 'soil, s=0.5', axs)
plot(resultB, 'r', 'peat, s=0.9', axs)
plot(resultB2, 'm', 'peat, s=0.5', axs)
plot(resultC, 'g', 'moss, s=0.9', axs)
plot(resultC2, 'lime', 'moss, s=0.5', axs)
finishPlot(axs)
plt.tight_layout()
plt.show()
In [ ]:
axs = getAxs()
resultAs = ewc_func(wrmA, poro=0.5, sl0=1.)
resultBs = ewc_func(wrmB, poro=0.7, sl0=1.)
resultCs = ewc_func(wrmC, poro=0.9, sl0=1.)
resultSurf = surf_func(1)
resultSurf2 = surf_func(.1)
plot(resultAs, 'b', 'soil', axs)
plot(resultBs, 'r', 'peat', axs)
plot(resultCs, 'g', 'moss', axs)
plot(resultSurf, 'goldenrod', 'surface 1.', axs)
plot(resultSurf2, 'k', 'surface 0.1', axs)
plt.show()
In [ ]:
In [ ]:
#wrmA2 = wrm_vangenuchten.VanGenuchten(m=0.1909, alpha=0.000545, sr=0.1, smoothing_interval=0.)
#wrmB2 = wrm_vangenuchten.VanGenuchten(m=0.4736, alpha=8.2e-6, sr=0.047058, smoothing_interval=1000.)
#wrmC2 = wrm_vangenuchten.VanGenuchten(m=0.275, alpha=0.00235, sr=0.055555, smoothing_interval=0.)
In [ ]:
kr = np.array([wrmA.k_relative(apc) for apc in pc])
kr2 = np.array([wrmA2.k_relative(apc) for apc in pc])
In [ ]:
s = np.array([wrmC.saturation(apc) for apc in pc])kr = np.array([wrmC.k_relative(apc) for apc in pc])
kr2 = np.array([wrmC2.k_relative(apc) for apc in pc])
plt.plot(s, kr,'b')
plt.plot(s, kr2,'r')
plt.xlim(0.9,1)
In [ ]:
s = np.array([wrmB.saturation(apc) for apc in pc])
kr = np.array([wrmB.k_relative(apc) for apc in pc])
kr2 = np.array([wrmB2.k_relative(apc) for apc in pc])
plt.plot(s, kr,'b')
plt.plot(s, kr2,'r')
plt.xlim(0.9999,1)
plt.ylim(0.9,1)
In [ ]:
ps1 = np.arange(0, 70000, 1000.)
ps2 = np.arange(70000, 110000, 100.)
ps3 = np.arange(110000, 150000, 1000.)
ps = np.concatenate((ps1, ps2, ps3))
def wc_func(wrm, poro, T, ps=ps):
# create the permafrost model
pm = permafrost_model_explicit_fpd.PermafrostModel(wrm)
wc = np.array([ewc.water_content(T,p,poro,pm) for p in ps])
sats = np.array([pm.saturations_Tp(T,p) for p in ps])
result = dict()
result['p'] = ps
result['wc'] = wc
result['T'] = T
result['porosity'] = poro
result['si'] = sats[:,2]
result['sl'] = sats[:,1]
result['sg'] = sats[:,0]
return result
def wc_getAxs(fig=None):
if fig is None:
fig = plt.figure(figsize=(14,4))
axs = []
axs.append(fig.add_subplot(121))
axs.append(fig.add_subplot(122))
axs[0].set_xlabel("pressure [Pa]")
axs[1].set_xlabel("pressure [Pa]")
axs[0].set_ylabel("water content [mol/m^3]")
axs[1].set_ylabel("gas saturation [-]")
return axs
def wc_plot(res, color, name, axs):
axs[0].plot(res['p'],res['wc'],'-x', color=color, label=name)
axs[0].plot(res['p'][1:],10000*(res['wc'][1:] - res['wc'][:-1])/(res['p'][1:] - res['p'][:-1]),'-x', color=color, label=name)
axs[1].plot(res['p'],res['sg'],'-x', color=color)
In [ ]:
# mineral soil
wrmA = wrm_vangenuchten.VanGenuchten(m=0.1909, alpha=0.000545, sr=0.1, smoothing_interval=100.)
resultA = wc_func(wrmA, 0.5, 273.65)
resultA2 = wc_func(wrmA, 0.5, 271.65)
axs = wc_getAxs()
wc_plot(resultA, 'b', 'soil, unfrozen', axs)
wc_plot(resultA2, 'c', 'soil, frozen', axs)
finishPlot(axs)
plt.tight_layout()
In [ ]:
ps1 = np.arange(-1000000, 70000, 5000.)
ps2 = np.arange(70000, 110000, 100.)
ps3 = np.arange(110000, 150000, 1000.)
ps = np.concatenate((ps1, ps2, ps3))
# peat soil
wrmB = wrm_vangenuchten.VanGenuchten(m=0.4736, alpha=8.2e-6, sr=0.047058, smoothing_interval=0)
resultB = wc_func(wrmB, 0.5, 273.65, ps)
resultB2 = wc_func(wrmB, 0.5, 271.65, ps)
axs = wc_getAxs()
wc_plot(resultB, 'r', 'peat, unfrozen', axs)
wc_plot(resultB2, 'm', 'peat, frozen', axs)
finishPlot(axs)
plt.tight_layout()
In [ ]:
# moss
wrmC = wrm_vangenuchten.VanGenuchten(m=0.275, alpha=0.00235, sr=0.055555, smoothing_interval=100.)
resultC = wc_func(wrmC, 0.5, 273.65)
resultC2 = wc_func(wrmC, 0.5, 271.65)
axs = wc_getAxs()
wc_plot(resultC, 'g', 'moss, unfrozen', axs)
wc_plot(resultC2, 'lime', 'moss, frozen', axs)
finishPlot(axs)
plt.tight_layout()
In [ ]:
wrm = wrmA
pm = permafrost_model_explicit_fpd.PermafrostModel(wrm)
p = 101647.904162
T = 273.120525
print pm.saturations_Tp(T,p)
In [ ]:
ps = np.linspace(101640, 101660, 10000)
wcs = np.array([ewc.water_content(T, ap, 0.5, pm) for ap in ps])
plt.plot(ps, wcs, '-x')
In [ ]:
Ts = np.linspace(T-1.e-4, T+1.e-4, 10000)
wcs = np.array([ewc.water_content(aT, p, 0.5, pm) for aT in Ts])
plt.plot(Ts, wcs)
In [ ]:
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.cm as cm
def plot_neighborhood(p,T,poro,wrm):
pm = permafrost_model_explicit_fpd.PermafrostModel(wrm)
ps = np.arange(p-1000.,p+1000.,5)
Ts = np.arange(T-1, T+1, 0.004)
es = np.zeros((len(ps),len(Ts)),'d')
wcs = np.zeros((len(ps),len(Ts)),'d')
for i in range(len(ps)):
for j in range(len(Ts)):
es[i,j],wcs[i,j] = ewc.ewc(Ts[j],ps[i],poro,pm)
PS,TS = np.meshgrid(ps,Ts)
centroid = ewc.ewc(T,p,poro,pm)
print p,T,centroid
fig = plt.figure()
ax1 = fig.add_subplot(121,projection='3d')
ax2 = fig.add_subplot(122,projection='3d')
ax1.plot_surface(PS.transpose(),TS.transpose(),es,cmap=cm.jet)
ax1.set_xlabel("pressure")
ax1.set_ylabel("temp")
ax1.set_zlabel("energy")
ax1.scatter([p],[T],[centroid[0]],s=100, c='white')
ax2.plot_surface(PS.transpose(),TS.transpose(),wcs, cmap=cm.jet)
ax2.set_xlabel("pressure")
ax2.set_ylabel("temp")
ax2.set_zlabel("wc")
ax2.scatter([p],[T],[centroid[1]],s=100, c='white')
In [ ]:
plot_neighborhood(85259.661,273.143,0.5,wrmA)
In [ ]:
plt.show()
In [ ]:
In [ ]: